/**************************************************************************
* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
*                                                                        *
* Author: The ALICE Off-line Project.                                    *
* Contributors are mentioned in the code where appropriate.              *
*                                                                        *
* Permission to use, copy, modify and distribute this software and its   *
* documentation strictly for non-commercial purposes is hereby granted   *
* without fee, provided that the above copyright notice appears in all   *
* copies and that both the copyright notice and this permission notice   *
* appear in the supporting documentation. The authors make no claims     *
* about the suitability of this software for any purpose. It is          *
* provided "as is" without express or implied warranty.                  *
**************************************************************************/
//
//  Xianguo Lu <lu@physi.uni-heidelberg.de>
//

#include <TFile.h>
#include <TH2I.h>
#include <TH2D.h>
#include <TList.h>
#include <THnSparse.h>
#include <TSystem.h>

#include "AliAnalysisManager.h"
#include "AliESDfriend.h"
#include "AliESDtrackCuts.h"
#include "AliESDInputHandler.h"
#include "AliTPCseed.h"
#include "AliTrackerBase.h"
#include "AliTRDseedV1.h"
#include "AliTRDtrackV1.h"
#include "AliCentrality.h"

#include "AliTRDdEdxUtils.h"
#include "AliTPCCosmicUtils.h"
#include "AliAnalysisTaskESDTRDsignal.h"

AliAnalysisTaskESDTRDsignal::AliAnalysisTaskESDTRDsignal()://const TString chunk):
  AliAnalysisTaskSE("ESDTRDsignal")
  //, fChunk(chunk)
  , fStream(0x0)
  , fESDEvent(0x0)

  , fHist(0x0)
{
  //
  //
  //

  DefineOutput(1, TTree::Class());

  fHist = new TList;
  fHist->SetOwner();//------>>> Very important!!!!!!!
  
  //fStream = new TTreeSRedirector(Form("TaskESDTRDsignalTree_%s.root", fChunk.Data()));
  fStream = new TTreeStream("tree");
  fTree=fStream->GetTree();
}

//____________________________________________________________
AliAnalysisTaskESDTRDsignal::~AliAnalysisTaskESDTRDsignal(){
  //
  // Destructor
  //
  //delete fStream;

  /* can not be deleted in chain, they are in input handler!
  delete fESDEvent;
  delete fESDfriend;
  delete fESDpid;
  */

  fHist->Clear();
  delete fHist;
}

void AliAnalysisTaskESDTRDsignal::UserExec(Option_t *)
{
  //
  //Execute
  // 
  Int_t mag = -999;
  if(!CutESD(mag))
    return;

  if(!CutEvent())
    return;

  const Int_t ntrk0 = fESDEvent->GetNumberOfTracks();
  AliESDtrack * trks[ntrk0];
  for(Int_t ii=0; ii<ntrk0; ii++){
    trks[ii] = 0x0;
  }
  Int_t nsel = 0;
  CutTrack(nsel, trks);
  if(!nsel)
    return;

  Double_t pbcen = -999;
  AliCentrality *esdCentrality = fESDEvent->GetCentrality();
  if(esdCentrality){
    pbcen = esdCentrality->GetCentralityPercentile("V0M");
  }
  
  Int_t run = fESDEvent->GetRunNumber();
  UInt_t timestamp = fESDEvent->GetTimeStamp();
  for(Int_t itrk=0; itrk<nsel ; itrk++){
    AliESDtrack * esdtrack = trks[itrk];

    if(!esdtrack->GetTRDntracklets())
      continue;
    
    Int_t ipar = -999;

    //in the tender it is the status bits of the tracks 14 - 16 for electrons, pions and protons
    Int_t tmpbit = 0;
    if(esdtrack->TestBit(BIT(14))) tmpbit+=1;
    if(esdtrack->TestBit(BIT(15))) tmpbit+=2;
    if(esdtrack->TestBit(BIT(16))) tmpbit+=4;

    if(tmpbit == 1) ipar = AliPID::kElectron; 
    if(tmpbit == 2) ipar = AliPID::kPion; 
    if(tmpbit == 4) ipar = AliPID::kProton;

    if(ipar<-990){
      continue;
    }
    
    const Double_t mass = AliPID::ParticleMass(ipar);

    Int_t charge = esdtrack->Charge();

    Double_t p0 = esdtrack->P();
    Double_t bg0 = p0/mass;
    Double_t ltbg0 = TMath::Log10(bg0);
    
    Double_t pTRD = -999;
    for(Int_t ilayer = 0; ilayer<6; ilayer++){
      const Double_t ptmp = esdtrack->GetTRDmomentum(ilayer);
      if(ptmp>1e-10){
        pTRD = ptmp;
        break;
      }
    }
    if(pTRD<0){
      continue;
    }

    Double_t bgTRD = pTRD/mass;
    Double_t ltbgTRD = TMath::Log10(bgTRD);

    Double_t tpcsig  = esdtrack->GetTPCsignal();
    Double_t trdsig = esdtrack->GetTRDsignal();
    Int_t trdnch  = esdtrack->GetTRDNchamberdEdx();
    Int_t trdncls = esdtrack->GetTRDNclusterdEdx();

    (*fStream)<<"tree"<<
      "run="<<run<<
      "timestamp="<<timestamp<<
      "pid="<<ipar<<
      "mag="<<mag<<
      "charge="<<charge<<
      "pbcen="<<pbcen<<

      "p0="<<p0<<
      "bg0="<<bg0<<
      "ltbgVtx="<<ltbg0<<

      "pTRD="<<pTRD<<
      "bgTRD="<<bgTRD<<
      "ltbg="<<ltbgTRD<<

      "esdtrack="<<esdtrack<<
      "tpcsig="<<tpcsig<<

      "trdsig="<<trdsig<<
      "trdnch="<<trdnch<<
      "trdncls="<<trdncls<<
      "\n";

  } 

  PostData(1, fTree);
}

void AliAnalysisTaskESDTRDsignal::Terminate(Option_t *) 
{
  //
  //do when Terminate
  //
}

//========================================================================================
//========================================================================================

Bool_t AliAnalysisTaskESDTRDsignal::CutEvent()
{
  Bool_t kpr0 = kTRUE;
  const AliESDVertex *vertex = fESDEvent->GetPrimaryVertexTracks();
  if(vertex->GetNContributors()<1) {
    // SPD vertex
    vertex = fESDEvent->GetPrimaryVertexSPD();
    if(vertex->GetNContributors()<1) {
      // NO GOOD VERTEX, SKIP EVENT 
      kpr0 = kFALSE;
    }
  }
  const Bool_t kpriv = kpr0 && ( fabs(fESDEvent->GetPrimaryVertex()->GetZ())<10 ); 

  return kpriv;

}

void AliAnalysisTaskESDTRDsignal::CutTrack(Int_t &nsel,  AliESDtrack *outtrk[])
{
  //
  //CutTrack
  //
  const Double_t etacut = 0.9;

  AliESDtrackCuts esdTrackCuts;//i.e. GetStandardITSTPCTrackCuts2010(kTRUE);

  //--->
  // TPC  
  esdTrackCuts.SetMinNCrossedRowsTPC(70);
  //esdTrackCuts.SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);

  esdTrackCuts.SetMaxChi2PerClusterTPC(4);
  esdTrackCuts.SetAcceptKinkDaughters(kFALSE);
  esdTrackCuts.SetRequireTPCRefit(kTRUE);
  // ITS
  esdTrackCuts.SetRequireITSRefit(kTRUE);
  /*
  esdTrackCuts.SetClusterRequirementITS(AliESDtrackCuts::kSPD,
                                        AliESDtrackCuts::kAny);
  // 7*(0.0026+0.0050/pt^1.01)
  esdTrackCuts.SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");

  esdTrackCuts.SetMaxDCAToVertexZ(2);
  esdTrackCuts.SetDCAToVertex2D(kFALSE);
  esdTrackCuts.SetRequireSigmaToVertex(kFALSE);
  */
  //---<
  
  //---
  esdTrackCuts.SetEtaRange(-etacut, etacut);
  //---------------------------------------------------------------<

  const TObjArray* seltrk = esdTrackCuts.GetAcceptedTracks(fESDEvent);
  nsel = seltrk->GetEntries();
  for(Int_t ii=0; ii<nsel; ii++){
    AliESDtrack *esdtrack = (AliESDtrack *)(seltrk->At(ii));
    outtrk[ii]=esdtrack;
  }
  delete seltrk;
}


Bool_t AliAnalysisTaskESDTRDsignal::CutESD(Int_t & mag)
{
  //
  //general protection
  //
  fESDEvent = 0x0;

  AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
  if(esdH){  
    fESDEvent = esdH->GetEvent();
  }

  if(!fESDEvent){
    printf("xlulog No ESD Event\n");
    exit(1);
  }
  
  //in unit of kg
  const Double_t dm = fESDEvent->GetMagneticField();
  if( fabs(fabs(dm)- 5) > 0.1 && fabs(fabs(dm)- 1) > 0.1){
    printf("xlulog strange Bfield! %f\n", dm);
    exit(1);
  }

  if(dm>0)
    mag = (int) (dm+0.5);
  else
    mag = (int) (dm-0.5);

  return kTRUE;
}

